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The phase behaviour and dynamics of molecular ionic liquids are studied using primitive models and exten- 
sive computer simulations. The models account for size disparity between cation and anion, charge location 
on the cation, and cation-shape anisotropy which are all prominent features of important materials such as 
room-temperature ionic liquids. The vapour-liquid phase diagrams are determined using high-precision Monte 
Carlo simulations, setting the scene for in-depth studies of ion dynamics in the liquid state. Molecular dynamics 
simulations are used to explore the structure, single-particle translational and rotational autocorrelation func- 
tions, cation orientational autocorrelations, self diffusion, viscosity and frequency-dependent conductivity The 
results reveal some of the molecular-scale mechanisms for charge transport, involving molecular translation, 
rotation, and association. 
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1. Introduction 

Molecular ionic liquids are low-melting point compounds made up of molecular cations and an- 
ions. Currently, the most important examples are room-temperature ionic liquids (RTILs). which 
have very low vapour pressures, solvate a wide range of polar and non-polar solutes, and have been 
shown to catalyse chemical transformations when used as reaction media [T . Typical cations (such 
as dialkylimidazolium) may exhibit both ionic and non-ionic characteristics. For instance, solutions 
of RTILs and water lead to the formation of micellar phases |2HB]. In terms of transport properties, 
RTILs are generally considered quite viscous and the diffusion constants are correspondingly low as 
compared to those in non-ionic molecular liquids. Experimentally, measurements of the frequency- 
dependent dielectric response yield insights on microscopic motions that lead to a change of polar- 
ization and, through the fundamental link between dielectric response and conductivity |7J[^, to 
the transport of charge. The dielectric spectra of many common RTILs have been measured j^HH]- 
Of particular note, Weingartner and colleagues found evidence for a contribution to the dielectric 
spectra of l-alkyl-3-methylimidazolium salts likely arising from cation rotations |15H17| . An analy- 
sis of dielectric spectra and nuclear-magnetic relaxation measurements highlighted motions on the 
order of tens of picoseconds (corresponding to real frequencies v ~ 10^^ Hz). Computer-simulation 
results have been used to resolve the different contributions to the dielectric spectrum, with the 
conclusion that the dielectric relaxation arising from molecular translations is faster than that from 
molecular rotations |15j . 

Atomistic simulations have yielded invaluable insights on the important molecular interactions 
and correlations in specific systems [18., 19j; there are far too many simulation studies to mention 
here, but a review by Maginn highlights the main achievements and outstanding challenges |20| . 
Low-frequency dynamics can be studied in atomistic simulations [21] but with some difficulty. 
The primary problem is that to calculate transport coefficients or frequency-dependent response 
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functions requires accurate calculations of appropriate time-correlation functions, but the rate 
of decay of these functions is characteristically long in RTILs and it is a challenge to calculate 
the long-time tails with acceptable signal-to-noise ratios [20 . Free-energy calculations and the 
construction of phase diagrams are other areas where the complexity of atomistic models may 
push the calculations out of reach. Therefore, to achieve a more comprehensive survey of dynamics 
over broad timescales - as well as the structure and phase behaviour - one might consider simplified 
models in which certain molecular characteristics can be tuned at will. 

Clearly, the main molecular characteristics of RTILs are that the ions are of different sizes, 
the charges are distributed unevenly over the ions (especially in typical cations), and the cations 
are normally far from spherical. Typical anions such as Cl~, BF4 , or PF,^, may be considered 
spherical in a first approximation, but cations such as dialkylimidazolium or dialkylpiperidinium are 
heterocyclic species. In the same way that a fluid of charged hard spheres, the restricted primitive 
model (RPM), serves as a basic model of simple molten salts, certain extensions to allow for the 
aforementioned molecular characteristics can be studied to gain insight on molecular ionic liquids. 
Malvaldi and Chiappe studied the dynamics in primitive models of RTILs consisting of dumbbell 
cations (each made up of two soft spheres, with one or both carrying a central charge) and simple 
soft-sphere anions |Z2] . It was shown that, even with such a simple model, certain experimental 
observations could be reproduced and rationalised in microscopic terms. For instance, in some 
imidazolium RTILs, the diffusion constant of the cation is greater than that of the anion |23l [21] ; 
atomistic simulations show that the cation diffuses preferentially by moving in the direction of 
the carbon in the 2 position on the five-membered ring (in between the nitrogens at the 1 and 3 
positions) pS]. This translational anisotropy is captured to some degree by Malvaldi and Chiappe 's 
dumbbell model where the charge is distributed equally between the two halves; at low temperature 
and high density, the cation diffuses faster than the anion. Spohr and Patey have performed a 
highly systematic and comprehensive survey of the microscopic structure and dynamics in models 
consisting of a large spherical cation carrying a point charge displaced from the centre, and a simple 
soft-sphere anion. By controlling the location of the cation charge and cation-anion size disparity, 
a number of experimental trends could be rationalised on the basis of ion-ion correlations |26fl28] . 
Specifically, many of the experimentally observed dependences of shear viscosity, conductivity, and 
diffusion constants on relative ion size, molecular charge distribution, and temperature could be 
mimicked with the simple models. Recent work on such models has even captured the effects of a 
polar solvent/impurity, such as water |29| . 

Regarding phase behaviour, the available thermodynamic data on RTILs show some interesting 
trends. For instance, the vapour-liquid critical temperature (Tc) - although difficult to measure 
directly in experiments because of molecular decomposition at high temperatures - appears to 
decrease with increasing molecular weight of the alkyl chains on the imidazolium cations, de- 
spite the growing van der Waals interactions favouring an increase in Tc jSOj. Martm-Betancourt 
et al. constructed coarse-grained models in which the cations are represented by charged hard 
spherocylinders, and performed high-precision Monte Carlo (MC) simulations to determine the 
coexistence envelopes and critical points [5T]. The critical temperature and density both decrease 
with increasing cation elongation, in qualitative correspondence with experimental data, due to the 
growing entropic role played by the steric bulk of the cation. Schroer and Vale surveyed fluid-fluid 
phase separations in solutions of imidazolium salts in a variety of polar and non-polar solvents, and 
showed that the phase diagrams could be analysed in terms of those for charged hard spheres with 
size asymmetry [32]. The main conclusion here was that the phase separation is mainly driven 
by Coulombic interactions, and that steric and dispersion interactions (which can be controlled 
systematically with substituents on the cations) modify the critical parameters. Models such as 
size-asymmetric charged hard spheres can be studied to yield systematic variations in critical pa- 
rameters with cation size; integral-equation studies point to some non-trivial effects P51 [51] . 

In this work, primitive models of molecular ionic liquids are constructed and their phase be- 
haviour and dynamical properties are examined. The models are chosen to reflect some of the key 
molecular characteristics, such as a size disparity between cation and anion, the location of the 
positive charge on the molecular cation, and the shape anisotropy of the cation. By constructing 
simple models possessing each of these characteristics, the essential relationships between molecu- 
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lar properties and dynamical behaviour can be identified. First, the vapour-Hquid phase diagrams 
are computed using high-precision MC simulations. Then, molecular dynamics (MD) simulations 
of the dense liquid are performed in order to evaluate the microscopic structure and the transport 
coefficients - diffusion constants, viscosity, and conductivity. The dependences of these properties 
on molecular architecture and temperature are briefly discussed. Finally, the frequency-dependent 
conductivity is considered and the contributions from cation translations and rotations are resolved 
by comparing the conductivity spectrum with the spectra of the translational and rotational- 
velocity autocorrelation functions; the orientational correlations of the cations are also considered. 
The possible relevance of these results to experimental measurements of the dielectric response is 
discussed. 

This paper is organised as follows. In section [2j the molecular models and parameters are 
defined. Monte Carlo simulations of phase behaviour are presented in section [3] focusing on the 
vapour-liquid coexistence envelope and the associated critical point. Molecular dynamics simu- 
lations of microscopic structure and (frequency-dependent) transport properties are discussed in 
section [4] Section [s] concludes the paper. 



2. Models 

The primitive models are presented in figure [l] In each case the anion is modelled as a single 
charged Lennard-Jones (LJ) sphere with energy parameter e, diameter a, charge —q, and mass 2m. 
In the symmetric dumbbell (SD) and asymmetric dumbbell (AD) models, each cation is modelled 
with two LJ spheres (with LJ parameters e and a, and mass m) fused together at a distance 
d = (t/\/2: in the SD model, each sphere carries a central charge of +q/2; in the AD model, 
only one sphere carries a central charge +q. In the corresponding symmetric-sphere (SS) and 
asymmetric-sphere (AS) models, each cation is formed from a single large LJ sphere with diameter 
(T_i_ , energy parameter e, and mass 2m. The diameter is chosen so that the excluded volume 
47rcri]_/3 is equal to that of the dumbbells defined above. The excluded volume of two dumbbells 
with sphere diameter a and separation d is given approximately by |35| 



TTCr 

~6~ 



8 + 12 



(2.1) 



With d = a/y/2, the effective cation diameter is ~ 1.3014(t. The charges on the cations are 
placed so that the shortest distance to the edge is always a/ 2: in the SS model, two charges +q/2 




Figure 1. Models of primitive liquids distinguished by the model for the cation: symmetric 
dumbbell (SD); asymmetric dumbbell (AD); symmetric sphere (SS); asymmetric sphere (AS); 
centred sphere (CS). In each case, the anion is modelled as a single charged LJ sphere with 
diameter a, charge —q, and mass 2m. The SD and AD cations are each modelled with two LJ 
spheres of diameter a and mass m, fused together at a distance d — a jsfi: in the SD model, each 
sphere carries a charge +q/2\ and in the AD model, only one sphere carries the charge ~\-q. The 
SS, AS, and CS cations are each modelled by a charged LJ sphere with diameter (t+ = 1.3014cr 
and mass 2m: the SS cation carries two charges +g/2 displaced symmetrically from the centre 
by 0.1507a"; the AS cation has a single charge -|-g located 0.1507ct from the centre; the CS 
cation has its charge +g at the centre. With these choices, the mass of each cation is 2m, the 
shortest distance from the positive charges(s) to the edge is a/2 (except for the CS model), and 
the excluded volume of each cation is the same. For the SD, AD, SS, and AS models, arrows 
indicate the unit vectors used to define the orientational autocorrelation function in section l473l 
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are displaced symmetrically from the centre by a distance (ct+ — a)/2 = 0.1507ct; in the AS model, 
one charge +g is placed 0.1507ct from the centre. Finally, in the centred-sphere (CS) model, the 
cation consists of a single charged LJ sphere with diameter (t+ ~ 1.3014ct, energy parameter e, 
mass 2m, and central charge +q. With this set of parameters, the excluded volumes for all cations 
are equal, the shortest distances from the charges to the edges are all equal to a/2 (except for the 
CS model), the charges on the cations and anions are ±g, and the cations and anions have equal 
masses of 2m. The interactions between ions are given by a sum of LJ and Coulombic (C) site-site 
interactions. The LJ interaction between two spheres i and j is given by 



/ \ 12 




(^) - 


(^) 







(2.2) 



where Uij = {ui + aj)/2, Tij is the centre-centre distance between spheres, and for simplicity the 
energy parameter e is the same for all interactions. The charge-charge interaction is 



(M3_ 

S Si'j 



(2.3) 



where qi is the charge on particle i, Sij is the distance between charges i and j, and £ = Atteq 
in which Eq is the dielectric permittivity of the vacuum. The reduced charge is defined as q* = 
\J q^ / Eae; in this work, there is given a fixed value of q* — 7. This is appropriate for monovalent 
ions (q — e) and realistic choices for the LJ energy parameter e/k^ ~ 600 K and the LJ diameter 
cr ~ 5 A. For each of the dumbbell models, the reduced moment of inertia of the cation is /* = 
I/ma'^ — 0.25; the same value is used for each of the spherical-cation models. Other reduced units 
are defined using e, a, and m: the reduced temperature is T* = k^T/e, where fee is Boltzmann's 
constant; the reduced ion-pair concentration p* = Na^ /V , where N is the number of ion pairs 
and V is the volume; the reduced pressure P* = Pa^/e; and the reduced time t* = t/r where 
T = ma^ / e . With the typical values of e and a given above, along with a characteristic salt 
molecular weight of 4m ~ 300 g mol~^, the basic unit of time is r ~ 1.94 ps. 

The phase diagrams of the model systems were determined using MC simulations; the simulation 
method and results are presented in section [3] Subsequently, the dynamical properties were deter- 
mined using MD simulations; the computational details and results are presented in section|4] In all 
cases, the systems were simulated in cubic boxes of side L and with periodic boundary conditions 
applied. The LJ interactions were truncated and shifted at L/2, and the Coulombic interactions 
were handled using the Ewald summation with conducting periodic boundary conditions |36| . 



3. Phase behaviour 

3.1. Monte Carlo simulations 

To identify the liquid region of the phase diagram and to focus the subsequent dynamical 
studies, vapour-liquid coexistence envelopes were determined using a MC technique. Wang-Landau 
simulations were performed in the grand-canonical ensemble according to the GCMCWL scheme 
described in reference |37]. Essentially, this is a flat-histogram sampling method that enables an 
iterative determination of the canonical partition function Q{N, V, T) as a function of N at fixed 
V and T. The conditions for phase coexistence are then easily determined using an equal-area 
criterion applied to the bimodal particle-number distribution p{N) oc Q{N,V,T), where z is 
the activity. Full details are given in reference |37| . The current simulations were carried out in 
cubic boxes with sides L = 11.75a and L = 13a for dumbbell and sphere models, respectively. The 
GCMCWL technique has been successfully applied to a number of 'tough' systems where a large 
degree of clustering is anticipated near the coexistence region, including charged soft spheres |37| 
and dipolar spheres |38| . 
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3.2. Vapour-liquid coexistence envelopes 

Vapour- liquid coexistence envelopes for the SD. AD. and AS systems are shown in figure [2] In 
order to assess the effects of the Coulombic interactions, results are also shown for an equimolar 
mixture of uncharged dumbbells and LJ spheres (the SD/AD model with g = 0). The critical 
parameters were estimated by fitting the simple scaling law 



P± = Pc 



At ± Bt 



(3.1) 



to near-critical coexistence densities in the vapour (p-) and liquid {p+) phases, where pc is the 
critical density, t = 1 — T/T^ , and /3 = 0.3265 is the 3D Ising order-parameter exponent. Ionic 
criticality is known to be Ising- like [35]. The fitted scaling laws are shown in figure |2] 
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Figure 2. (Color on-line) Phase diagrams of selected primitive models of ionic fluids: (circles) 
SD model; (squares) AD model; (diamonds) AS model; (triangles) SD/AD model with g = 0. 
Open symbols indicate sub-critical coexistence points, filled symbols indicate the critical points. 



and the lines show the fits from the scaling law in equation (3.1 1. The dashed line indicates the 



temperature range along the p* = 0.27 isochore studied using molecular dynamics simulations. 

The critical parameters are given in table [T] By comparing the charged and uncharged cases, 
it is clear that the Coulombic interactions result in a large increase in the critical temperature 
and a substantial decrease in the critical density. A low critical density is one characteristic of 
an ionic critical point [40J, and with the current parameters, the Coulombic interactions are very 
strong. To compare the current results with those for the simple case of charged hard spheres - 
the restricted primitive model (RPM) - table [T| also shows the packing fraction 77 at the critical 
point. For the systems studied here, the volumes of the ions are estimated by approximating them 
as hard particles. For the SD and AD models, the volume of a dumbbell is 



v{d) 
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(3.2) 



The key point is that rj^ — 0.04 for all charged models, including the RPM. The comparison 
of critical temperatures is not that straightforward. The RPM critical temperature is k-eTc = 
0.05069(2)g^/f a |39], where —q^ /8a is the Coulombic interaction between a cation and an anion 
at contact. For the primitive models studied here, the minimum in the interaction potential between 
two oppositely charged LJ spheres is rmi„ ~ 0.9556cr and Wmin = — 49.63e. The reduced critical 
temperatures are all in the region of T* — 6 and so the corresponding 'ionic' temperature is 
fceT'c/lumml — 0.12. This is considerably higher than the RPM value, refiecting the soft cores of 
the particles and the additional LJ attractions. Model AD has the highest critical temperature, 
refiecting the fact that there are a large number of cation-anion arrangements where the charge 



33602-5 



G.C. Ganzenmuller, P.J. Camp 



TabiB 1 . Vapour-liquid critical parameters for selected primitive models of ionic fluids (including 
the RPM [39]) and for an uncharged system of dumbbells and LJ spheres (SD/AD model with 
q = 0). T* = ksTc/e is the reduced critical temperature, p* — pcfJ^ is the reduced critical ion- 
pair density, and ?7c is the critical packing fraction calculated from p* assuming hard particles 
of dimension a. The figures in brackets are the estimated uncertainties in the final digits. 
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SD 


5.91(2) 


0.0255(1) 


0.0385(3) 


AD 


6.22(2) 


0.0281(2) 


0.0424(3) 


AS 


5.80(2) 


0.0247(1) 


0.0414(2) 


SD/AD (q = 0) 


1.970(8) 


0.0671(6) 


0.1013(9) 


RPM [31] 


0.05069(2)(j*2 


0.0395(1) 


0.0414(1) 



separation is about la and the Coulombic interaction is strong. The number of such arrangements 
is reduced when the cation charge is split between two centres, as in the SD model, because now 
the anion has to be in the right place to satisfy two site-site interactions. For model AS, only 
one arrangement of the anion and cation gives a charge-charge separation of about la, and so the 
critical temperature is relatively low. 

Knowledge of the vapour-liquid coexistence region facilitates the choice of the temperature 
and density for MD simulations in the liquid phase. The aim is to reach low temperatures (and 
hence low vapour pressures) while at the same time staying within the liquid region of the phase 
diagram. The choice p* = 0.27 gives access to the liquid region of the phase diagram. Moreover, 
for a typical molecular ionic liquid with cr = 5 A and ion mass 2m ~ 150 g mol^^, this density 
equates to a mass density of about 1 g cm^'^, which closely corresponds to experimental values. 
The temperature range 1.5 ^ T* ^ 3.0 will be considered, as the upper limit is well within the 
liquid region for all models, while the lower limit approaches the region of spinodal decomposition. 



4. Liquid-phase structure and dynamics 

4.1. Molecular dynamics simulations 

MD simulations were performed using the LAMMPS package [TO I42| . Systems of iV = 125 ion 
pairs were simulated in a cubic box with edge length L — 7.736cr giving a reduced ion-pair density 



of p* = 0.27, as advertised in section 3.2 The dynamical equations of motion were integrated 
using the velocity- Verlet scheme with a timestep of 6t* — 0.0025. Simulations were started from a 
CsCl lattice structure and equilibrated at the desired temperature using simple velocity rescaling 
over 10^ timesteps. Production runs were then performed in the canonical (NVT) ensemble over 
7.2 X 10^ timesteps (equivalent to around 35 ns) with weak coupling to a Berendsen thermostat |36| . 
Ion positions and velocities, and all components of the stress tensor, were saved at intervals of 4 



Tabis 2. Reduced pressure P* along the isochore p* — 0.27 as a function of temperature, for 
each of the charged models. 
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timesteps for post-processing. The production runs were divided into 12 blocks, and statistical 
uncertainties were estimated by assuming block averages to be statistically independent. Table [2] 
shows the values of the pressure along the isochore p* — 0.27 for each of the charged models. At 
temperatures T* < 2.5, the pressure in at least one of the systems is negative, indicating that these 
state points may well be mechanically unstable and/or within a binodal region. Figure [2] confirms 
that these low temperatures are approaching the vapour-liquid coexistence region. 

4.2. Structure 

The microscopic structure in the liquids was examined using the radial distribution functions 
(RDFs) gafsir) {a, /3 = +, — ) [33]. These functions are given by 



V 



I Nff \ 



(4.1) 



where = N- = N, and r^j is the distance between the centres of mass of particles i and j. MD 
simulation results at T* — 2.5 are shown in figure [3j this temperature was selected in order to be 
sure of simulating a single dense liquid phase. Firstly, the cation-cation RDF (5++) and anion-anion 




Figure 3. (Color on-line) Pair-correlation functions in the liquid phase at p* = 0.27 and T* — 

2.5: (a) cation-cation RDF, (ji++(r); (b) cation-anion RDF, g-i (r); (c) anion-anion RDF, g (r). 

In all cases, the separation r is that between centres of mass. 

RDF {g ) show primary peaks at around r ~ 1.5cr for all models. The cation-anion RDF (g_| ) 

shows strong correlations at much shorter distances, but there is a marked difference between the 
spherical models (SS, AS, CS) and the dumbbell models (SD, AD), reflecting the shapes of the 
cations. The spherical models show one primary peak at r ~ l.lcr, roughly equal to the distance of 
closest approach (cr + (t+)/2 ~ 1.15cr. The dumbbells models exhibit a short-range double-peaked 
feature: the peak near r ~ Icr is clearly due to the anion sitting in the 'valley' between the two 
spheres of the cation dumbbell; the peak near r ~ 1.3(T roughly corresponds to the anion being 
located near the dumbbell axis, at a distance a + d/2 ~ 1.21cr from the cation centre of mass. This 
picture is in good correspondence with work on similar models by Malvaldi and Chiappe |22| . 



4.3. Dynamics 

A range of dynamical properties has been calculated. The frequency-dependent conductivity, 
k{uj), is related to the charge current deflned by 



(4.2) 



where the sum runs over all charge sites (of which there are Nq), Avi(t) is the velocity of charge 
i in its molecular centre-of-mass frame, and v™'"(t) is the corresponding molecular centre-of-mass 
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velocity. (Obviously, Av = for the anions.) The frequency-dependent conductivity is given by [JS] 

oo 

Kiuj) = + iK"{L0) = j {3{t) ■ J(0)) cxp [-iujt)dt, (4.3) 



where i = \/— 1. and k'{uj) and k"(u!) are the real (in-phase) and imaginary (out-of-phase) parts, 
respectively. The static conductivity is k = k{0). The contributions of different single-particle 
motions to the conductivity spectrum can be determined by calculating ion translational velocity 
autocorrelation functions Cy (t) ^43) , and a cation intramolecular velocity autocorrelation function 
CAit): 

Cyit) = (vr^(i).vr°(0)), (4.4) 
CA{t) = (Av,(t).Av,(0)). (4.5) 

These two functions capture the dynamics of the molecular centre of mass [Ct,(t)] and the in- 
tramolecular rotations [CA(i)l • The ion diffusion constants are obtained from the standard 
Green-Kubo relationship 

C30 

D= \ I Cy(t)dt. (4.6) 



3 



The reorientational dynamics of the cations are examined by calculating the correlation function 

Ce(t) = (e,(t)-e,(0)), (4.7) 

where e{t) is the orientation unit vector along the cylindrical symmetry axis of the cation. An 
associated correlation time can be defined by 

oo 

Te = j Ce{t)dt. (4.8) 


Finally, the viscosity of the fluid is calculated using the Green-Kubo formula involving the auto- 
correlation function of the off-diagonal elements of the stress tensor Iij.y{t): 

oo 

V = TTTT^ I (n,j,(t)n,j,(0))di. (4.9) 



VkBT 







The static transport properties are considered first, namely the diffusion constants (£)), shear 
viscosity (77), and static conductivity (k). In reduced units, these properties are defined by D* = 
Dt / , rf = r]a^ /tT, and k* — kt/£. These properties were measured for each charged model 
as functions of temperature. In accord with Eyring's activated-dynamics picture (44!, it turns out 
that the data are fitted adequately with Arrhenius relations of the form 

i?cxexp(-SaABT), r;cxexp(^a/fcBT), k cx exp (-i;a/A:BT), (4.10) 

where iJa is an activation energy associated with that particular transport property. The logarithms 
of the transport properties are plotted as functions of 1/T* in figure |4] along with Arrhenius-law 
fits. Simulations at T* = 1.5 gave erratic and anomalous results, probably due to the system being 
within either the liquid-solid or vapour-liquid coexistence region; these data points are omitted from 
the plots. Arrhenius-like behaviour is observed over the temperature range considered. In general, 
for a given temperature, the cation and anion diffusion constants obey Z?+ < Z3_ due to the larger 
steric bulk of the cations. The dumbbells models (SD and AD) exhibit higher values of -D+ than 
the sphere models (SS, AS, CS). This may be attributed to the elongated shape of the dumbbell 
cations, facilitating the motion along their long axes. For the dumbbell cations, a symmetric charge 
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Figure 4. (Color on-line) Arrhenius plots of transport properties at p* — 0.27 and 1.75 ^ T* ^ 
3.00: (a) cation diffusion constant D+; (b) anion diffusion constant D_; (c) shear viscosity 77; 
(d) static conductivity k. Results are plotted in reduced units defined in the text. 



distribution gives a higher value of _D_|_ than an asymmetric charge distribution, which probably 
arises from there being less strong interactions (on average) between an SD cation and its nearest- 
neighbour anions. For the spherical cations, the reverse is true, although the magnitude of the 
effect is smaller. In general, £)_ mirrors the trends seen in , indicating strong association 
between oppositely charged ions. Using the molecular parameters given in section [2] the reduced 
values D* ~ 0.03-0.10 correspond to real values D ~ 0.3-1.0 x 10~^ m^ s"-'^, which are a little 
higher than typical results for RTILs. The shear viscosity t] shows the reverse trend to D± , with 
the dumbbells models being least viscous; this can be rationalised qualitatively on the basis of a 
Stokes-Einstein relationship of the form D (x I/77 The reduced values ij* ~ 3-5 equate to real 
values of 77 ~ 3-5 x 10^^ Pa s. These are much smaller than the viscosities of real RTILs, which 
are normally on the scale of tens of mPa s. The trends in the static conductivity mirror those in 
D± , which can be rationalised with a generalised Nernst-Einstein relationship of the form [13] 

'^=l^{D+ + D_)il~A), (4.11) 

where A characterises the deviation from pure Nernst-Einstein behaviour, which is observed when 
there are no cross correlations between cation and anion motions. Such correlations lead to a 
reduction in k and hence A > 0. For the data shown in figure|4] A ~ 0.75, a large value that arises 
from very strong cation-anion correlations. The reduced values k* ~ 0.1-0.2 correspond to real 
values of K ~ 6-12 S m~^, which are somewhat higher than those seen in experiments on RTILs. 
In summary, then, for dumbbell cations, a symmetric charge distribution gives higher diffusion 
constants, lower shear viscosity, and higher static conductivity. The dependences of transport 
properties with spherical cations are less sensitive to the charge distribution. 

The activation energies associated with the static transport properties are reported in table |3] 
In general, the dumbbell-cation models exhibit lower activation energies. This correlates with the 
picture presented above, with the elongated models being more capable of breaking out of transient 
molecular cages by moving along their long axes. When converted to real energy units - using the 
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Tabl6 3. Activation energies obtained from Arrhenius-law fits to the temperature dependence 
of the cation diffusion constant (D+), the anion diffusion constant the shear viscosity {rj), 

the static conductivity (k), and an orientational correlation time (re). Figures in brackets are 
estimated uncertainties in the final digits. 





Activation energies £'a/e 


Model 


D+ 






K 




SD 


3.28(13) 


3.71(13) 


1.27(21) 


0.823(80) 


2.452(43) 


AD 


3.63(12) 


3.636(96) 


1.47(10) 


1.532(93) 


2.938(51) 


SS 


4.51(14) 


4.580(89) 


1.806(73) 


1.70(14) 


1.301(39) 


AS 


3.98(26) 


4.10(25) 


1.55(14) 


1.47(11) 


1.748(37) 


CS 


4.79(18) 


4.546(38) 


1.725(62) 


1.51(12) 





molecular parameters given in section [2]- the activation energies for diffusion are in the region of 
20 kJ mol~^. which is a typical value for some RTILs [13 . The activation energies for viscosity 
and conductivity are in the region of 5-10 kJ moP^, which are the right order of magnitude but a 
little too low. Experimental results for RTILs show significant deviations from Arrhenius behaviour 
and are usually fitted with a Vogel-Fulcher-Tammann equation of the form Aexp [B/{T — Tq)], so 
these are not directly comparable to the Arrhenius-law fits presented here. Nonetheless, B/k-B is 
typically of order 500-1000 K, which in the current reduced units equates to around B/e ~ 1-2; 
this is at least the same order of magnitude as the activation energies presented in table [3) 

The dynamics are now considered in more detail for a single state point, p* = 0.27 and T* = 
2.5. This temperature is chosen such that each of the models is in its liquid phase, away from 
the coexistence region. Figure |5] shows the conductivity spectra for all five models. The 

spectra of the SD and AD models are characterised by broad shoulders in the real parts, k'{oj), at 
frequencies ~ 10-20, and single peaks in the imaginary parts, k"{uj), at a frequency w* ~ 30. 
The spectrum for the CS model is particularly simple, with the real part exhibiting a peak that 
is characteristic of ionic liquids ^45n47] . The low-frequency, negative portion of the imaginary part 
signals that ion translations may contribute to the low-frequency complex dielectric spectrum 



0.30 I — I I I I I I 1 1 — I — I — I I I I I 1 1 — I — I I I I I I r 




Figure 5. Real and imaginary parts of the conductivity spectrum, Hi [u)) (solid lines) and ^"(a;) 
(dashed lines), respectively, for primitive models at p* — 0.27 and T* = 2.5: (a) SD model; (b) 
AD model; (c) SS model; (d) AS model; (e) CS model. 
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e(a;) = 1 + 47riK(cj)/w [5T]. The spectrum of the SS model looks very similar to that of the CS 
model, showing that a symmetrical displacement of charge on the cation does not strongly affect the 
conductivity. This is natural, because these two models share the same centre of mass and the same 
centre of charge; the displacement of charge should only affect the degree of interactions between 
cations and anions, and is not expected to lead to new qualitative features in k{uj). Very interesting 
features arise in the spectrum of the AS model. Here, the real part contains two, very-well resolved 
peaks at reduced frequencies of around uj* = 5 and uj* = 20. 

To help identify the molecular motions responsible for the features in figure [6] shows the 

real part along with the spectra of Cv{t) (for each of the ions) and that of CA{t): 

oo 

C{u}) = 2 j C{t) cos {ujt)dt. (4.12) 



Recall that Cy{t) describes the centre-of-mass motions of the ions, while intramolecular rotations 
of the cations are described by CA{t)- The spectra are basically densities of states for vibrations 
and librations within the solvation shell. Of course, these are single-particle properties, and in the 
following discussion, the cross-correlations between cations and anions (that lead to deviations 
from the Nernst-Einstein law) are ignored. Nonetheless, the single-particle properties will yield 
some insight on the collective dynamics that lead to the conductivity spectrum. For the purposes 
of comparison, each spectrum is normalised to unit area. For the SD and AD models, the cation 
translations, anion translations, and cation rotations occur on similar timescales. As a result, 
k'{uj) exhibits a broad peak over the relevant range of frequencies. For the SS and AS models, the 
rotational spectra are much less broad than those for the SD and AD models. For the AS model, 
the peak in the rotational spectrum coincides with the low-frequency peak in k'{uj), indicating that 
this feature of the conductivity arises from cation rotations. Of course, this feature is absent from 
k'{uj) for the SS model, because the rotation of a symmetrical distribution of positive ions leads to 
no net transport of charge. 

The peaks in the rotational spectra in figure ^ show that the spherical cations (SS and AS) 




Figure 6. (Color on-line) Real part of the conductivity spectrum [^'(a;)] compared with the spec- 
tra of the centre-of-mass velocity autocorrelation functions for the cations and anions [C„(aj)], 
and the rotational velocity autocorrelation function for the cations [CAC^i^)], for primitive models 
at p* = 0.27 and T* = 2.5: (a) SD model; (b) AD model; (c) SS model; (d) AS model. In each 
case, the plots are normalised so that the frequency integral for each spectrum is equal to unity. 



33602-1 1 



G.C. Ganzenmuller, P.J. Camp 



librate slower than the dumbbell cations, despite their having the same excluded volumes and equal 
moments of inertia. This difference is attributed to the so-called charge arm, described by Kobrak 
and Sandalow [48 . If a single charge q is displaced by a distance Iq from the cation centre of 
mass, then the librational frequency of rotation is proportional to Iq/I [4S3; Iq is called the charge 
arm. Recall that in the present models, the charges are set a certain distance from the edge of 
the cations, and so they may be at different distances from the centres of mass. For the dumbbell 
models Iq — a/2y/2 ~ 0.3536(7, while for the spherical models Iq = 0.1507cr. Keeping the moment 
inertia /* ~ 0.25 for the AD model and equating the ratio Iq/I for the AD and AS models gives 
/* = 0.1066 for the AS model. Conductivity spectra and the spectra of the velocity autocorrelation 
functions have been calculated for the SS and AS models with /* = 0.1066 (data not shown) and 
the rotational peak does indeed shift to higher frequency, coinciding with the peak frequencies for 
the SD and AD models. Correspondingly, the low-frequency peak in k'(cj) shifts up in frequency 
and merges with the centre-of-mass peak. 

The single-particle motions detailed here occur in the region of a;* ~ 5-20, which equate to real 
frequencies v ~ 0.4-1.6 THz. These are quite close to the natural timescales in RTILs: cation ro- 
tations in l-alkyl-3-methylimidazolium salts occur on the timescale of tens of picoseconds |15ffT7] . 
Shim and Kim present k{u)) from computer simulations of l-ethyl-3-methylimidazolium hexaflu- 
orophosphate that exhibits a main peak at around 15 THz and a low-frequency shoulder below 
about 10 THz [2T|; it is tempting to speculate that the low-frequency feature arises from cation 
rotations. 




t* l/T* 

Figure 7. (Color on-line) (a) Orientational autocorrelation function Ce{t) for primitive models 
at p* = 0.27 and T* — 2.5. (b) Arrhenius plot of the corresponding decay times for 1.75 ^ 
T* ^ 3.00. Results are plotted in reduced units defined in the text. 

To explore cation reorientations further, figure [t] (a) shows Ce{t) for the dumbbell-cation and 
spherical-cation models, at p* = 0.27 and T* = 2.5. Firstly, the dumbbell models exhibit a 
monotonous decay in Ce{t), while the spherical models show a negative portion. The spherical 
models may undergo a rapid reversal in orientation on the timescale of t ~ 0.6r due to the absence 
of steric hindrance. Secondly, the correlations die away more slowly with dumbbells than with 
spheres, due to the packing effects that restrict the reorientation of the anisotropic dumbbells. 
Finally, an asymmetric charge distribution, as in the AD and AS models, leads to a slower decay 
of Ce{t) as compared to the SD and SS models, respectively. This may be due to the anisotropy of 
the cation-anion interactions when a single, whole charge is localised near the edge of the cation, 
rather than being spread out over the molecule. To summarise, anisotropies in the short-range re- 
pulsions and electrostatic interactions lead to longer orientational correlation times. The decay time 
T* — Te/r has been calculated as a function of temperature for each of the models; the results are 
presented on Arrhenius plots in figure [t] (b) . The associated activation energies from the Arrhenius 
law Te oc exp (£^a/fcB2^) are reported in table [3j The activation energies for the dumbbell models 
are greater than those for the spherical models, once again reflecting the steric barriers to rotation 
of dumbbells arising from their anisotropic shapes. An asymmetric charge distribution gives rise to 
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marginally higher barriers than does a symmetric one, reflecting the effect of anisotropy of cation- 
anion interactions. The decay times span the range r* ~ 0.1-0.8, which in real units correspond 
to the picosecond timescale; these are realistic values, at least for simple molecular liquids. When 
converted to real units, the activation energies are of order 10 kJ moP^, which as mentioned above, 
are at least of the right order of magnitude for molecular ionic liquids, but a bit too low. 



5. Conclusions 

Simple models of molecular ionic liquids were constructed in order to study the effects of 
ion-size disparity, cation-charge distribution, and cation-shape anisotropy. The cations were made 
up of charged soft spheres, with attractive interactions included: one type of cation consisted of 
a single sphere carrying either one charge at the centre or off-centre, or two charges displaced 
symmetrically from the centre; the other type of cation was a dumbbell, with either one or two 
of the constituent spheres carrying a charge. The anions were charged soft spheres. Such simple 
models admit in-depth studies of phase behaviour and dynamics in simulations of reasonable length, 
using normal computational facilities. The vapour-liquid coexistence curves were determined using 
Wang-Landau simulations and the resulting critical parameters were compared with those of the 
restricted primitive model, the simplest representation of an ionic fluid. The presence of attractive 
interactions increases the critical temperature above that expected from electrostatic interactions 
alone, while the critical density remains roughly constant and characteristically low. 

Dynamics simulations were then conducted along a liquid-state isochore appropriate to typical 
ambient conditions. The microscopic structure has been studied using pair correlation functions, 
and rationalised in terms of the probable packing of the ions, and the resulting electrostatic inter- 
actions. The temperature dependences of the ion diffusion constants, shear viscosity, conductivity, 
and orientational correlation time are Arrhenius-like, at least over the temperature range con- 
sidered here. The magnitudes of these properties, and the associated activation energies, when 
converted to real units, are not too dissimilar from those measured experimentally, although in 
general the diffusion constants and conductivity are too high, and the viscosity is too low. This is to 
be expected of such simple models with comparatively little 'molecular roughness'. Coarse graining 
usually results in an effective decrease of molecular-scale friction, implying enhanced diffusion and 
conductivity, and reduced viscosity. This is a general drawback of coarse-grained models [49| which 
cannot always be remedied just by tuning the effective interactions; for instance, molecular-scale 
noise and friction can be restored with integration schemes similar to those used in dissipative 
particle dynamics simulations |50.. The variations in the transport coefficients between different 
models were rationalised in terms of microscopic correlations. The dynamics at a particular liquid- 
phase state point (high density, low temperature) were investigated in more detail by looking at 
the frequency-dependent conductivity, and the spectra of centre-of-mass and rotational velocity au- 
tocorrelation functions. For the models with dumbbell cations, the centre-of-mass and rotational 
motions occur on similar timescales, although only the centre-of-mass motions contribute to the 
conductivity spectrum if the positive charge is distributed symmetrically over the molecule. For 
the models with spherical cations and off-centre charges, the centre-of-mass and rotational motions 
occur on different timescales which, with an asymmetrical distribution of charge on the cation, lead 
to well-resolved features in the conductivity spectrum. The decay of orientational correlations was 
rationalised with reference to anisotropics of short-range repulsions and electrostatic interactions. 
Despite the simplicity of the models studied here, it is possible that some of the features observed 
in the detailed dynamics may also be observable in experimental studies of molecular ionic liquids. 
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OaaoBa noBefliHKa ra ^MHaiviiKa b npuMirnBHux MOfle/iax 

MO/ieKy/15ipHO-iOHHI/1X pi^i^H 

r.C. FaHseHMKD/i/ie 

iHCTMxyT fli/iHaMiKM uj BUflKOTpM Bail nx npoL4eciB TOBapi/icTBa 0payHro4>epa (iHCTHxyT EpHcra Maxa), 
OpaCiGypr, HiMeHHWHa 

LLlKO/ia xImIT, yHlBepci/iTeT EfliH5ypra, EfllHSypr, Bejii/iKo6pMTaHiiR 

floc.nifl>KeHO cfjaBOBy noBefllHKy I flMHaMiKy Mo;ieKy;i5ipHO-ioHHHX pifli/iH 3 flonoMororo npi/iMiri/iBHi/ix MOfle- 
f\e\A Ta MaciuTa6Hwx KOMn'raxepHnx MOfle/iKDBaHb. BuKopucraHi MOflejii BpaxoBytoxb pisHi/iLiK) b poBMipax 
KaxioHa xa ahioHa, no/io>KeHHiR sapjRfly Ha KaxioHax i aHlBoxponlKi cfjopMH KaxioHa, mo e BM3HaHajibHi/iMH 
B/iacxHBOcxaMM ioHHMX piflMH Hpi/i KlMHaxHi/ix xeM nepaxy pax. BucoKoxoHHe MOflejiHDBaHHa MexoflOM MoHxe- 
Kapjio Bi/iKopi/icxaHO fl/ia no5yflOBH cJds^obi/ix fliarpaw piflUHa-ras, jikI b nofla/ibiuoMy BUKopucxoByraxbca 

fl/ia BHBHeHHJI fll/lHaMiKM ioHHI/lX piflHH. 3a flOnOMOfOKD MOfle/lKDBaHHiR MOJlEKyjia pHOtO flMHaMlKOKD flOC/lifl>Ke- 

Ho cxpyKxypy, xpaHC/iJiLiiiiHi xa o5epxa/ibHi aBXOKopejiiRL4iMHi 4>yHKi4iT, KaxioHHi opieHxaL4iMHi aBXOKopejiiRL4iV, 

CaMOfll/l4)y3iK), B'aBKicXb Xa HaCXOXHy 3ajie>KHicXb npOBlflHOCxi. Oxpi/IMaHi peSyjlbXaXH Bl/lJIB/iaKDXb HI/13Ky MO- 

/leKy/iiqpHwx MexaHiBMlB nepeHocy BapiqfliB, BKJiKDHaKDHH MOJieKyjiapHy xpaHCJinL4iKD, o5epxaHH5i xa acoL4iai4iKi. 
K/iiOHOBi c/iosa: ionni pifli/iHi/i, nepexift pifli/iHa-ra3, fti/inaMiKa, KOMn'ioTepHe MOfte/iioBaHHft 
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